clear; close all; clc;

format long

% Note: the variables v1, v2, v3 are used to determine which regions have
% solid lines versus dashed lines (decreasing vs increasing parts of the
% schedule). They were adjusted manually.

%% Benchmark case

cases = 1;
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure 8: Interest rate correspondence for Spain (b = 15%)
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4350;v2=4375;v3=4595;
ig=2;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0.8],'LineWidth',2.00,'LineStyle','--'); hold on
ig=2;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
ig=2;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0.8],'LineWidth',2.00,'LineStyle','--');
ig=2;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',10)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',10)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',10)
xlim([5 35])
ylim([-1 12])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,'low growth','high growth');
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.5]; % [width height]
paper_position = [0 0 5.0 2.5]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rcorresp.pdf';
saveas(gcf,filename)

% Figure 9.b : Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - Inherited debt b
bb1 = 0.15;
bb2 = 0.17;
ib1v  = find(bvec<bb1,1,'last');
ib2v  = find(bvec<bb2,1,'last');
figure
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib2v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib2v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib2v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib2v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([-2 32])
ylim([-1 13])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$b$=',num2str(round(100*bb2,0)),'\%'], ...
             ['$b$=',num2str(round(100*bb1,0)),'\%']);
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.0]; % [width height]
paper_position = [0 0 5.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rsch_b.pdf';
saveas(gcf,filename)

% Figure 11.a: Policy functions and equilibrium spreads for Spain - Debt issuance (n), low growth
figure
ig=1;is=2;iz=floor(Nz/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p1=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50); hold on
ig=1;is=1;iz=floor(Nz/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p2=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 20])
ylim([5 45])
hleglines = [p1(1) p2(1)];
leg = legend(hleglines,'good sunspot','bad sunspot');
set(leg,'Interpreter','LaTex','Fontsize',8,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_npol_gL.pdf';
saveas(gcf,filename)

% Figure 11.b: Policy functions and equilibrium spreads for Spain - Debt issuance (n), high growth
figure
ig=2;is=2;iz=floor(Nz/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p3=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50); hold on
ig=2;is=1;iz=floor(Nz/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p4=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 20])
ylim([5 45])
hleglines = [p3(1) p4(1)];
leg = legend(hleglines,'good sunspot','bad sunspot');
set(leg,'Interpreter','LaTex','Fontsize',8,'Location','SouthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_npol_gH.pdf';
saveas(gcf,filename)

% Figure 11.c: Policy functions and equilibrium spreads for Spain - Spreads (ρ − r∗), low growth
figure
ig=1;is=2;iz=floor(Nz/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p1=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50); hold on
ig=1;is=1;iz=floor(Nz/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p2=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9)
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^{*}$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 20])
ylim([0 5])
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rpol_gL.pdf';
saveas(gcf,filename)

% Figure 11.d: Policy functions and equilibrium spreads for Spain - Spreads (ρ − r∗), high growth
figure
ig=2;is=2;iz=floor(Nz/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p3=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50); hold on
ig=2;is=1;iz=floor(Nz/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p4=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9)
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^{*}$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 20])
ylim([0 5])
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rpol_gH.pdf';
saveas(gcf,filename)

%% Varying the probability of the bad sunspot: psB

cases = [1 2];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure 9.a : Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - The role of expectation pB
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=3528;v2=3550;v3=3795;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
%title(['Interest rate correspondence, b = ',num2str(round(100*bb,0)),'\%'],'Interpreter','LaTex','Fontsize',23)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([-2 32])
ylim([-1 13])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$p_B$=',num2str(round(param(12,2),2))], ...
             ['$p_B$=',num2str(round(param(12,1),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.0]; % [width height]
paper_position = [0 0 5.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_pBs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_pBs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table 4: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Case 1: pB=',num2str(round(param(12,1),2)),'  (benchmark) \n'])
fprintf(['Case 2: pB=',num2str(round(param(12,2),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Case 1               \t Case 2                \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying maturity (delta)

cases = [3 1 4];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);
ibv= zeros(ncases,1);
DD = 1.00;

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);    
    load(['OUTPUT_',num2str(cases(ic)),'/bvec.txt'])     
    ibv(ic)  = find(bvec<DD*param(2,ic),1,'last');             
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure 7.c : Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - Maturity δ (b/δ = 100%)
figure
ig=1;is=1; xx = 100*nsch(:,ibv(1),ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>20;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<20)|(yy>35);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ibv(1),ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0.5 0],'LineWidth',3.00,'LineStyle',':'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ibv(1),ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0.5 0],'LineWidth',1.0,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ibv(1),ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0.5 0],'LineWidth',3.00,'LineStyle',':');
ig=1;is=1;p4=plot(100*nsch(v3:end,ibv(1),ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0.5 0],'LineWidth',1.0,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ibv(2),ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ibv(2),ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ibv(2),ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ibv(2),ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ibv(2),ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ibv(3),ig,is,3); yy = 100*(Rsch(:,ig,is,3)-(param(1,3)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p9=plot(100*nsch(1:v1,ibv(3),ig,is,3),100*(Rsch(1:v1,ig,is,3)-(param(1,3)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p10=plot(100*nsch(v1:v2,ibv(3),ig,is,3),100*(Rsch(v1:v2,ig,is,3)-(param(1,3)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p11=plot(100*nsch(v2:v3,ibv(3),ig,is,3),100*(Rsch(v2:v3,ig,is,3)-(param(1,3)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p12=plot(100*nsch(v3:end,ibv(3),ig,is,3),100*(Rsch(v3:end,ig,is,3)-(param(1,3)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([-2 32])
ylim([-1 13])
hleglines = [p1(1) p5(1) p9(1)];
leg = legend(hleglines,...
             ['$\delta$=',num2str(round(param(2,1),2)),' (' ,num2str(round(1/param(2,1),1)),' years)'], ...
             ['$\delta$=',num2str(round(param(2,2),2)),' (' ,num2str(round(1/param(2,2),1)),' years)'], ...
             ['$\delta$=',num2str(round(param(2,3),2)),' (' ,num2str(round(1/param(2,3),1)),' years)']); 
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','SouthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.0]; % [width height]
paper_position = [0 0 5.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_deltas_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_deltas.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.3: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (2): delta=',num2str(round(param(2,1),2)),'\n'])
fprintf(['Column (1): delta=',num2str(round(param(2,2),2)),'  (benchmark) \n'])
fprintf(['Column (3): delta=',num2str(round(param(2,1),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (2)           \t Column (1)           \t Column (3)            \n')
fprintf('-------------------                \t -------------------- \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying persistence of low-growth state (pL)

cases = [1 5];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);              
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.8a: Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - Probability of remaining in low growth (pL)
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=4305;v2=4330;v3=4572;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([5 40])
ylim([-1 13])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$p_L$=',num2str(round(param(6,1),2))], ...
             ['$p_L$=',num2str(round(param(6,2),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_pgLs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_pgLs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.3: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (1): pgLs=',num2str(round(param(6,1),2)),'  (benchmark) \n'])
fprintf(['Column (4): pgLs=',num2str(round(param(6,2),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (1)           \t Column (4)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the recovery rate (kappa)

cases = [1 6];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);    
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure 9.7b: Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - Recovery rate (κ)
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=4200;v2=4227;v3=4452;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
%title(['Interest rate correspondence, b = ',num2str(round(100*bb,0)),'\%'],'Interpreter','LaTex','Fontsize',23)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([5 40])
ylim([-1 13])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$\kappa$=',num2str(round(param(5,1),2))], ...
             ['$\kappa$=',num2str(round(param(5,2),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_kappas_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_kappas.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
         
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.3: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (1): kappa=',num2str(round(param(5,1),2)),'  (benchmark) \n'])
fprintf(['Column (5): kappa=',num2str(round(param(5,2),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (1)           \t Column (5)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the dispersion of idiosyncratic shocks (sigma)

cases = [1 7];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.8c: Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - St. deviation of transitory shock (σ)
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=4347;v2=4378;v3=4610;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',3.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.0,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',3.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.0,'LineStyle',':');
v1=4580;v2=4628;v3=4870;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([5 40])
ylim([-1 13])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$\sigma$=',num2str(round(100*param(11,1),1)),'\%'], ...
             ['$\sigma$=',num2str(round(100*param(11,2),1)),'\%']); 
set(leg,'Interpreter','LaTex','Fontsize',8,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_sigmas_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_sigmas.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.3: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (1): sigma=',num2str(round(param(11,1),3)),'  (benchmark) \n'])
fprintf(['Column (6): sigma=',num2str(round(param(11,2),3)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (1)           \t Column (6)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the default cost in the high-growth state (ybarH)

cases = [1 8];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.8d: Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - Cost of default in high growth (ϕ(gH))
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=4258;v2=4286;v3=4463;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
%title(['Interest rate correspondence, b = ',num2str(round(100*bb,0)),'\%'],'Interpreter','LaTex','Fontsize',23)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([5 40])
ylim([-1 13])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$\phi (g_{H})$=',num2str(round(param(4,1),3))], ...
             ['$\phi (g_{H})$=',num2str(round(param(4,2),3))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_ybarHs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_ybarHs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.3: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (1): ybarH=',num2str(round(param(4,1),3)),'  (benchmark) \n'])
fprintf(['Column (7): ybarH=',num2str(round(param(4,2),3)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (1)           \t Column (7)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the persistence of the high-growth stat (pH)

cases = [1 9];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                 
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.8e: Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - Probability of remaining in high growth (pH)
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=4840;v2=4870;v3=5128;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
%title(['Interest rate correspondence, b = ',num2str(round(100*bb,0)),'\%'],'Interpreter','LaTex','Fontsize',23)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([5 40])
ylim([-1 13])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$p_H$=',num2str(round(param(7,1),2))], ...
             ['$p_H$=',num2str(round(param(7,2),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','SouthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_pgHs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_pgHs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.3: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (1): pH=',num2str(round(param(7,1),2)),'  (benchmark) \n'])
fprintf(['Column (8): pH=',num2str(round(param(7,2),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (1)           \t Column (2)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the discount factor (beta)

cases = [1 10];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                 
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.8f: Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - Discount factor (β)
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=5040;v2=5070;v3=5367;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
%title(['Interest rate correspondence, b = ',num2str(round(100*bb,0)),'\%'],'Interpreter','LaTex','Fontsize',23)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([5 40])
ylim([-1 13])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$\beta$=',num2str(round(param(8,1),2))], ...
             ['$\beta$=',num2str(round(param(8,2),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','SouthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_betas_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_betas.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.3: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (1): beta=',num2str(round(param(8,1),2)),'  (benchmark) \n'])
fprintf(['Column (9): beta=',num2str(round(param(8,2),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (1)           \t Column (9)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the difference between high and low growth rates (gH-gL)

cases = [1 11];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.8g: Interest rate spreads for Spain in the low growth state
% (b = 15%): Comparative statics - Difference between growth rates (gH − gL)
bb = 0.15;
ib1v  = find(bvec<bb,1,'last');
figure
v1=4965;v2=4990;v3=5178;
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
v1=4347;v2=4378;v3=4607;
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
%title(['Interest rate correspondence, b = ',num2str(round(100*bb,0)),'\%'],'Interpreter','LaTex','Fontsize',23)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([5 40])
ylim([-1 13])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$\Delta^{g}$=',num2str(round(100*(param(10,1)-param(9,1)),0)),'\%'], ...
             ['$\Delta^{g}$=',num2str(round(100*(param(10,2)-param(9,2)),0)),'\%']); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_gHgLs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_gHgLs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.3: Simulation moments: Spain \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (1):  gH-gL=',num2str(round(param(10,1)-param(9,1),2)),'  (benchmark) \n'])
fprintf(['Column (10): gH-gL=',num2str(round(param(10,2)-param(9,2),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (1)           \t Column (10)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the probability of the bad sunspot: psB
% ALTERNATIVE CALIBRATION (DEBT = 70% of GDP)

cases = [12 13];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.3: Interest rate correspondence for southern Europe (low-growth state)
bb1 = 0.12;
bb2 = 0.11;
ib1v  = find(bvec<bb1,1,'last');
ib2v  = find(bvec<bb2,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>5;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<5)|(yy>12);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib2v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>5;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<5)|(yy>12);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib2v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib2v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib2v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib2v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',10)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',10)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',10)
xlim([4 16])
ylim([0 11])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$b$=',num2str(round(100*bb2,1)),'\%'], ...
             ['$b$=',num2str(round(100*bb1,1)),'\%']);
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.5]; % [width height]
paper_position = [0 0 5.0 2.5]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_Rsch_b_alt_calib.pdf';
saveas(gcf,filename)

% Figure F.4a: Policy functions and equilibrium spreads for Spain - Debt issuance (n), low growth
figure
ig=1;is=2;iz=floor(Nz/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p1=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50); hold on
ig=1;is=1;iz=floor(Nz/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p2=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 15])
ylim([5 45])
hleglines = [p1(1) p2(1)];
leg = legend(hleglines,'good sunspot','bad sunspot');
set(leg,'Interpreter','LaTex','Fontsize',8,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_npol_gL_alt_calib.pdf';
saveas(gcf,filename)

% Figure F.4b: Policy functions and equilibrium spreads for Spain - Debt issuance (n), high growth
figure
ig=2;is=2;iz=floor(Nz/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p3=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50); hold on
ig=2;is=1;iz=floor(Nz/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p4=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 15])
ylim([5 45])
hleglines = [p3(1) p4(1)];
leg = legend(hleglines,'good sunspot','bad sunspot');
set(leg,'Interpreter','LaTex','Fontsize',8,'Location','SouthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_npol_gH_alt_calib.pdf';
saveas(gcf,filename)

% Figure F.4c: Policy functions and equilibrium spreads for Spain - Spreads (ρ − r∗), low growth
figure
ig=1;is=2;iz=floor(Nz/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p1=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50); hold on
ig=1;is=1;iz=floor(Nz/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p2=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9)
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^{*}$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 15])
ylim([0 5])
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rpol_gL_alt_calib.pdf';
saveas(gcf,filename)

% Figure F.4d: Policy functions and equilibrium spreads for Spain - Spreads (ρ − r∗), high growth
figure
ig=2;is=2;iz=floor(Nz/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p3=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50); hold on
ig=2;is=1;iz=floor(Nz/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p4=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9)
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^{*}$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 15])
ylim([0 5])
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rpol_gH_alt_calib.pdf';
saveas(gcf,filename)

filename = 'Results/Table_pBs_alt_calib.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.1: Simulation moments: Spain low-debt calibration \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Case 1: pB=',num2str(round(param(12,1),2)),'  (benchmark) \n'])
fprintf(['Case 2: pB=',num2str(round(param(12,2),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Case 1               \t Case 2                \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off